First we import our tidied salamander morphology & trait dataset.
# Import tidied salamander eye size and trait data
salamanders <- data.frame(read.csv("../Data/Tidy/salamanders_tidy.csv", header=TRUE, na.strings=c("", "NA", " ","<1")))
# Quick look at data structure
str(salamanders)
## 'data.frame': 424 obs. of 39 variables:
## $ X : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Genus_Species : chr "Ambystoma_amblycephalum" "Ambystoma_amblycephalum" "Ambystoma_annulatum" "Ambystoma_annulatum" ...
## $ Order : chr "Caudata" "Caudata" "Caudata" "Caudata" ...
## $ Suborder : chr "Salamandroidea" "Salamandroidea" "Salamandroidea" "Salamandroidea" ...
## $ Family : chr "Ambystomatidae" "Ambystomatidae" "Ambystomatidae" "Ambystomatidae" ...
## $ Subfamily : chr NA NA NA NA ...
## $ Genus : chr "Ambystoma" "Ambystoma" "Ambystoma" "Ambystoma" ...
## $ Species : chr "amblycephalum" "amblycephalum" "annulatum" "annulatum" ...
## $ SVL_mm : int 73 60 80 80 88 60 56 102 103 115 ...
## $ Mass_g : num 11.58 9.34 8.97 8.39 12.12 ...
## $ rootmass : num 2.26 2.11 2.08 2.03 2.3 ...
## $ ED_right_mm : num 4.11 3.82 3.1 3.45 3.4 2.98 2.24 5.25 4.77 4.95 ...
## $ ED_left_mm : num 3.96 3.63 3.17 3.51 3.6 2.58 2.17 5.2 4.83 4.87 ...
## $ eyemean : num 4.04 3.73 3.13 3.48 3.5 ...
## $ CD_right_mm : num 2.93 2.42 2.06 2.34 2.2 1.7 1.52 3.37 3.56 3.62 ...
## $ CD_left_mm : num 2.94 2.11 1.99 2.44 2.43 1.7 1.75 3.32 3.33 4.11 ...
## $ cormean : num 2.94 2.27 2.02 2.39 2.31 ...
## $ Gill_Presence : chr "yes" "yes" "no" "no" ...
## $ Aquatic : chr "yes" "yes" "no" "no" ...
## $ Semiaquatic : chr "no" "no" "no" "no" ...
## $ Terrestrial : chr "no" "no" "no" "no" ...
## $ Scansorial : chr "no" "no" "no" "no" ...
## $ Fossorial : chr "no" "no" "yes" "yes" ...
## $ Subfossorial : chr "no" "no" "no" "no" ...
## $ Cave_dwelling : chr "no" "no" "no" "no" ...
## $ Nocturnal : chr NA NA "yes" "yes" ...
## $ Diurnal : chr NA NA "no" "no" ...
## $ Paedomorphic : chr "yes" "yes" "no" "no" ...
## $ Metamorphosizing : chr "no" "no" "yes" "yes" ...
## $ Direct_Development: chr "no" "no" "no" "no" ...
## $ Free_living : chr "yes" "yes" "yes" "yes" ...
## $ No_free_living : chr "no" "no" "no" "no" ...
## $ Lentic_water : chr "yes" "yes" "yes" "yes" ...
## $ Lotic_water : chr "no" "no" "no" "no" ...
## $ Both_water : chr "no" "no" "no" "no" ...
## $ No_larva : chr "no" "no" "no" "no" ...
## $ Females_larger : chr "yes" "yes" "yes" "yes" ...
## $ Males_larger : chr "no" "no" "no" "no" ...
## $ No_SSD : chr "no" "no" "no" "no" ...
Next, we import the amphibian tree from Jetz and Pyron (2019) that has been pruned to match the salamander species in this dataset.
#Import pruned phylogeny
caudatatree <- read.nexus("../Data/Tidy/caudata-tree.nex")
#Plot tree
plot.phylo(caudatatree, show.tip.label=TRUE, cex = 0.7)
First, we can take a preliminary look at all of our data so that we can see how it looks and whether we notice any obvious outliers that we might want to go back and check out. Here we plot out our measurements from the right side vs. the left side of each specimen so that we can make sure that we have consistent measurements throughout the dataset. We use the package plotly to create interactive plots - hover over a point to see the name of the species measured.
We expect that measurements of left and right eyes and corneas should be similar within specimens. We can plot out our data to test for this and make sure we don’t see any glaring issues. We can also fit a regression, which we would expect to have a slope near 1 and an intercept near 0 if left and right eye measurments are similar.
####right eye diameter vs. left eye diameter ####
#linear model (ordinary least squares model)
RLeye.lm <- lm(ED_right_mm ~ ED_left_mm, data = salamanders)
summary(RLeye.lm)
##
## Call:
## lm(formula = ED_right_mm ~ ED_left_mm, data = salamanders)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.53420 -0.13491 -0.01135 0.12950 0.67196
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.015122 0.035243 0.429 0.668
## ED_left_mm 1.010149 0.009162 110.253 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2045 on 383 degrees of freedom
## (39 observations deleted due to missingness)
## Multiple R-squared: 0.9695, Adjusted R-squared: 0.9694
## F-statistic: 1.216e+04 on 1 and 383 DF, p-value: < 2.2e-16
#plot
RLeye.plot <- ggplot(salamanders, aes(x = ED_left_mm, y = ED_right_mm, text = Genus_Species)) +
geom_point(alpha = 0.9) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
geom_abline(slope = coef(RLeye.lm)[[2]], intercept = coef(RLeye.lm)[[1]]) #plots OLS fit
#interactive plot
ggplotly(RLeye.plot)
As we would expect, measurements of the left and right eye are highly correlated, with a slope close to 1 and an intercept close to 0. We can do the same test for the corneal measurements.
####right cornea diameter vs. left cornea diameter ####
#linear model (ordinary least squares model)
RLcornea.lm <- lm(CD_right_mm ~ CD_left_mm, data = salamanders)
summary(RLcornea.lm)
##
## Call:
## lm(formula = CD_right_mm ~ CD_left_mm, data = salamanders)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.71444 -0.11655 -0.00936 0.13772 0.78095
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.10579 0.03248 3.257 0.00122 **
## CD_left_mm 0.97409 0.01177 82.742 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2133 on 422 degrees of freedom
## Multiple R-squared: 0.9419, Adjusted R-squared: 0.9418
## F-statistic: 6846 on 1 and 422 DF, p-value: < 2.2e-16
#plot
RLcor.plot <- ggplot(salamanders, aes(x = CD_left_mm, y = CD_right_mm, text = Genus_Species)) +
geom_point(alpha = 0.9) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
geom_abline(slope = coef(RLcornea.lm)[[2]], intercept = coef(RLcornea.lm)[[1]]) #plots OLS fit
#interactive plot
ggplotly(RLcor.plot)
Again, we find that left and right corneal measurements are highly correlated, with a slope close to 1 and an intercept close to 0. Together, these results indicate that measurements were consistent between the left and right eyes of a specimen, which is good! For looking at interspecific allometry, we will calculate species means from these data so that we can use phylogenetic comparative methods.
Here, we will use phylogenetic least-squares (PGLS) regressions to see how eyes are scaling with body size across different species of salamanders with different ecologies. PGLS analyses are similar to normal linear regressions, but also account for the non-independence of data collected from close relatives due to shared evolutionary history.
In these analyses, we will log eye size and body size variables, as this is standard for allometric relationships.
First, we will calculate the mean measurements across specimens for each species. We need to make sure in doing this that we omit rows with missing data for a variable of interest. In looking at the data, we can see that we have many more measurements for cornea diameter than for eye diameter
#find species with missing eye measurements
eye.missing <- salamanders %>%
filter(is.na(eyemean)) %>%
droplevels() %>%
group_by(Family,Genus_Species) %>%
summarise(n = n())
#create RMarkdown table of species missing eye data
kable(eye.missing, caption = "Specimens without eye size data") %>%
kable_styling(full_width = F) %>%
collapse_rows(columns = 1, valign = "top") %>%
scroll_box(height = "500px")
| Family | Genus_Species | n |
|---|---|---|
| Amphiumidae | Amphiuma_means | 3 |
| Amphiumidae | Amphiuma_tridactylum | 3 |
| Cryptobranchidae | Andrias_davidianus | 2 |
| Cryptobranchidae | Andrias_japonicus | 3 |
| Cryptobranchidae | Cryptobranchus_alleganiensis | 4 |
| Plethodontidae | Eurycea_nana | 3 |
| Plethodontidae | Eurycea_rathbuni | 3 |
| Proteidae | Necturus_alabamensis | 1 |
| Proteidae | Necturus_beyeri | 2 |
| Proteidae | Necturus_lewisi | 3 |
| Proteidae | Necturus_maculosus | 3 |
| Salamandridae | Pseudobranchus_striatus | 3 |
| Sirenidae | Siren_intermedia | 3 |
| Sirenidae | Siren_lacertina | 3 |
#find species missing cornea measurements
cor.missing <- salamanders %>%
filter(is.na(cormean)) %>%
droplevels() %>%
group_by(Family,Genus_Species) %>%
summarise(n = n())
#species missing cornea measurements
print(cor.missing)
## # A tibble: 0 × 3
## # Groups: Family [0]
## # … with 3 variables: Family <chr>, Genus_Species <chr>, n <int>
## # ℹ Use `colnames()` to see all variable names
We then find species means for various morphological parameters. There are some specimens that lack data for eye size but have data for svl, mass, and corneal diameter (these measurements are present for all specimens). However, when eye size measurements are missing, they are consistently missing for all specimens measured within a species. This means we don’t have to worry about missing data for one specimen corrupting the mean across several specimens for eye size (it is either present for all specimens in a species or absent for all specimens in a species). So we can make just one dataframe of species averages, and some species will just lack eye size data.
# species means for cornea data
av.sal <- salamanders %>%
mutate_if(is.character, as.factor) %>%
group_by(Genus_Species) %>%
summarise(eye_av = mean(eyemean),
cor_av = mean(cormean),
svl_av = mean(SVL_mm),
rootmass_av = mean(rootmass),
n = n())
## Merge data means with other column info, keeping only trait columns
av.sal <- merge(av.sal, salamanders[match(unique(salamanders$Genus_Species), salamanders$Genus_Species), ], by="Genus_Species", all.x = TRUE, all.y = FALSE) %>%
select(Genus_Species, eye_av, cor_av, svl_av, rootmass_av, n, Order, Suborder, Family, Subfamily, Genus, Species, Aquatic, Semiaquatic, Terrestrial, Scansorial, Fossorial, Subfossorial, Cave_dwelling, Nocturnal, Diurnal, Paedomorphic, Metamorphosizing, Direct_Development, Free_living, No_free_living, Lentic_water, Lotic_water, Both_water, No_larva, Females_larger, Males_larger, No_SSD)
In order to fit a PGLS regression in caper, we first need to make a comparative data object that includes our dataset and our phylogeny. Note that the gilled Pleurodeles waltl will be dropped from this, as there can only be one row for each species.
#Make row names of dataset the species names (so it will match phylogeny tips)
rownames(av.sal) <- av.sal$Genus_Species
#check that names match in dataframe and tree
name.check(phy = caudatatree, data = av.sal)
## $tree_not_data
## character(0)
##
## $data_not_tree
## [1] "Pleurodeles_waltl_gilled"
#use caper function to combine phylogeny and data into one object (this function also matches species names in tree and dataset)
sal.comp <- comparative.data(phy = caudatatree, data = av.sal,
names.col = Genus_Species, vcv = TRUE,
na.omit = FALSE, warn.dropped = TRUE)
#check for dropped tips or dropped species
sal.comp$dropped$tips #phylogeny
## character(0)
sal.comp$dropped$unmatched.rows #dataset
## [1] "Pleurodeles_waltl_gilled"
#PGLS eye diameter vs. cube root of mass
pgls_edmass <- pgls(log10(eye_av) ~ log10(rootmass_av),
data = sal.comp,
lambda = "ML", #uses Maximum Likelihood estimate of lambda
param.CI = 0.95)
#diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_edmass)
par(mfrow = c(1, 1))
#Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval
lambda.edmass <- pgls.profile(pgls_edmass, "lambda")
plot(lambda.edmass)
#print model output
summary(pgls_edmass)
##
## Call:
## pgls(formula = log10(eye_av) ~ log10(rootmass_av), data = sal.comp,
## lambda = "ML", param.CI = 0.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0180704 -0.0036954 -0.0000669 0.0035026 0.0228160
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.845
## lower bound : 0.000, p = 3.0531e-11
## upper bound : 1.000, p = 2.8592e-06
## 95.0% CI : (0.658, 0.942)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.398691 0.042603 9.3584 2.22e-16 ***
## log10(rootmass_av) 0.773017 0.041458 18.6458 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.006815 on 136 degrees of freedom
## (14 observations deleted due to missingness)
## Multiple R-squared: 0.7188, Adjusted R-squared: 0.7167
## F-statistic: 347.7 on 1 and 136 DF, p-value: < 2.2e-16
#plot
edmass.plot <- ggplot(data = av.sal, aes(x = log10(rootmass_av), y = log10(eye_av), text = Genus_Species)) +
geom_point(alpha = 0.7) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
geom_abline(slope = coef(pgls_edmass)[[2]], intercept = coef(pgls_edmass)[[1]]) #plots OLS fit
#interactive plot
ggplotly(edmass.plot)
#PGLS eye diameter vs. SVL
pgls_edsvl <- pgls(log10(eye_av) ~ log10(svl_av),
data = sal.comp,
lambda = "ML", #uses Maximum Likelihood estimate of lambda
param.CI = 0.95)
#diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_edsvl)
par(mfrow = c(1, 1))
#Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval
lambda.edsvl <- pgls.profile(pgls_edsvl, "lambda")
plot(lambda.edsvl)
#print model output
summary(pgls_edsvl)
##
## Call:
## pgls(formula = log10(eye_av) ~ log10(svl_av), data = sal.comp,
## lambda = "ML", param.CI = 0.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0173372 -0.0021347 0.0006661 0.0041418 0.0154655
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.749
## lower bound : 0.000, p = 3.9968e-14
## upper bound : 1.000, p = 1.1559e-08
## 95.0% CI : (0.506, 0.898)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.839890 0.079508 -10.564 < 2.2e-16 ***
## log10(svl_av) 0.782916 0.039652 19.745 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.005834 on 136 degrees of freedom
## (14 observations deleted due to missingness)
## Multiple R-squared: 0.7414, Adjusted R-squared: 0.7395
## F-statistic: 389.9 on 1 and 136 DF, p-value: < 2.2e-16
#plot
edsvl.plot <- ggplot(data = av.sal, aes(x = log10(svl_av), y = log10(eye_av), text = Genus_Species)) +
geom_point(alpha = 0.7) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
geom_abline(slope = coef(pgls_edsvl)[[2]], intercept = coef(pgls_edsvl)[[1]])
#interactive plot
ggplotly(edsvl.plot)
#PGLS eye diameter vs. cornea diameter
pgls_edcd <- pgls(log10(eye_av) ~ log10(cor_av),
data = sal.comp,
lambda = "ML", #uses Maximum Likelihood estimate of lambda
param.CI = 0.95)
#diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_edcd)
par(mfrow = c(1, 1))
#Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval
lambda.edcd <- pgls.profile(pgls_edcd, "lambda")
plot(lambda.edcd)
#print model output
summary(pgls_edcd)
##
## Call:
## pgls(formula = log10(eye_av) ~ log10(cor_av), data = sal.comp,
## lambda = "ML", param.CI = 0.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0069593 -0.0012901 -0.0000271 0.0017229 0.0066368
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.347
## lower bound : 0.000, p = 6.9782e-05
## upper bound : 1.000, p = < 2.22e-16
## 95.0% CI : (0.096, 0.648)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.169381 0.012943 13.087 < 2.2e-16 ***
## log10(cor_av) 0.962429 0.021166 45.471 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.002278 on 136 degrees of freedom
## (14 observations deleted due to missingness)
## Multiple R-squared: 0.9383, Adjusted R-squared: 0.9378
## F-statistic: 2068 on 1 and 136 DF, p-value: < 2.2e-16
#plot
edcd.plot <- ggplot(data = av.sal, aes(x = log10(cor_av), y = log10(eye_av), text = Genus_Species)) +
geom_point(alpha = 0.7) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
geom_abline(slope = coef(pgls_edcd)[[2]], intercept = coef(pgls_edcd)[[1]])
#interactive plot
ggplotly(edcd.plot)
#PGLS cornea diameter vs. cube root of mass
pgls_cdmass <- pgls(log10(cor_av) ~ log10(rootmass_av),
data = sal.comp,
lambda = "ML", #uses Maximum Likelihood estimate of lambda
param.CI = 0.95)
#diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_cdmass)
par(mfrow = c(1, 1))
#Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval
lambda.cdmass <- pgls.profile(pgls_cdmass, "lambda")
plot(lambda.cdmass)
#print model output
summary(pgls_cdmass)
##
## Call:
## pgls(formula = log10(cor_av) ~ log10(rootmass_av), data = sal.comp,
## lambda = "ML", param.CI = 0.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.031038 -0.005368 -0.000095 0.004781 0.041013
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.832
## lower bound : 0.000, p = 8.4377e-15
## upper bound : 1.000, p = 5.6595e-09
## 95.0% CI : (0.683, 0.922)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.184679 0.057335 3.2211 0.001566 **
## log10(rootmass_av) 0.586098 0.051971 11.2773 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01006 on 150 degrees of freedom
## Multiple R-squared: 0.4588, Adjusted R-squared: 0.4552
## F-statistic: 127.2 on 1 and 150 DF, p-value: < 2.2e-16
#plot
cdmass.plot <- ggplot(data = av.sal, aes(x = log10(rootmass_av), y = log10(cor_av), text = Genus_Species)) +
geom_point(alpha = 0.7) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
geom_abline(slope = coef(pgls_cdmass)[[2]], intercept = coef(pgls_cdmass)[[1]])
#interactive plot
ggplotly(cdmass.plot)
#PGLS cornea diameter vs. SVL
pgls_cdsvl <- pgls(log10(cor_av) ~ log10(svl_av),
data = sal.comp,
lambda = "ML", #uses Maximum Likelihood estimate of lambda
param.CI = 0.95)
#diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_cdsvl)
par(mfrow = c(1, 1))
#Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval
lambda.cdsvl <- pgls.profile(pgls_cdsvl, "lambda")
plot(lambda.cdsvl)
#print model output
summary(pgls_cdsvl)
##
## Call:
## pgls(formula = log10(cor_av) ~ log10(svl_av), data = sal.comp,
## lambda = "ML", param.CI = 0.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.028530 -0.006068 0.001092 0.006296 0.043279
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.872
## lower bound : 0.000, p = < 2.22e-16
## upper bound : 1.000, p = 1.4035e-07
## 95.0% CI : (0.752, 0.942)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.664239 0.120396 -5.5171 1.475e-07 ***
## log10(svl_av) 0.537185 0.051402 10.4507 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01098 on 150 degrees of freedom
## Multiple R-squared: 0.4213, Adjusted R-squared: 0.4175
## F-statistic: 109.2 on 1 and 150 DF, p-value: < 2.2e-16
#plot
cdsvl.plot <- ggplot(data = av.sal, aes(x = log10(svl_av), y = log10(cor_av), text = Genus_Species)) +
geom_point(alpha = 0.7) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
geom_abline(slope = coef(pgls_cdsvl)[[2]], intercept = coef(pgls_cdsvl)[[1]])
#interactive plot
ggplotly(cdsvl.plot)
#PGLS svl vs. cube root of mass
pgls_svlmass <- pgls(log10(svl_av) ~ log10(rootmass_av),
data = sal.comp,
lambda = "ML", #uses Maximum Likelihood estimate of lambda
param.CI = 0.95)
#diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_svlmass)
par(mfrow = c(1, 1))
#Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval
lambda.svlmass <- pgls.profile(pgls_svlmass, "lambda")
plot(lambda.svlmass)
#print model output
summary(pgls_svlmass)
##
## Call:
## pgls(formula = log10(svl_av) ~ log10(rootmass_av), data = sal.comp,
## lambda = "ML", param.CI = 0.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0190565 -0.0039221 0.0004095 0.0038784 0.0182303
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.928
## lower bound : 0.000, p = < 2.22e-16
## upper bound : 1.000, p = 7.3054e-07
## 95.0% CI : (0.855, 0.970)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.617161 0.033287 48.582 < 2.2e-16 ***
## log10(rootmass_av) 1.004421 0.025367 39.596 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.005716 on 150 degrees of freedom
## Multiple R-squared: 0.9127, Adjusted R-squared: 0.9121
## F-statistic: 1568 on 1 and 150 DF, p-value: < 2.2e-16
#plot
svlmass.plot <- ggplot(data = av.sal, aes(x = log10(rootmass_av), y = log10(svl_av), text = Genus_Species)) +
geom_point(alpha = 0.7) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
geom_abline(slope = coef(pgls_svlmass)[[2]], intercept = coef(pgls_svlmass)[[1]])
#interactive plot
ggplotly(svlmass.plot)
Here, we plot absolute eye size (cornea diameter) and relative eye size (the residuals of cornea diameter vs. cube root of mass) onto the salamander phylogeny. We will use the ggtree package for plotting. Again, the gilled Pleurodes waltl will not be included in these plots (as it was not included in the PGLS regressions).
First we add PGLS residuals to the dataframe and resort the data to match the species order of the phylogeny.
#Add residuals from PGLS fits to dataframe
#extract pgls residuals
pglsres.cdmass <- as.data.frame(residuals(pgls_cdmass))
#name residuals
colnames(pglsres.cdmass) <- "pglsres.cdmass"
#add column for genus_species
pglsres.cdmass$Genus_Species <- rownames(pglsres.cdmass)
#merge residuals with original data by rowname
sal.res <- inner_join(av.sal, pglsres.cdmass, by = "Genus_Species")
#export dataset for use later
write.csv(sal.res, file = "../Data/Tidy/salamanders_residuals.csv")
#make dataframe for labeling phylogeny figure
sal.phy <- sal.res %>%
#make column labeling tip labels that match phylogeny
mutate(tip = Genus_Species) %>%
#put genus and species in separate columns
separate(Genus_Species, c("genus", "species"), sep = "_", extra = "drop") %>%
#add tip labels with regular text
mutate(labels = as.factor(paste(genus, species, sep = " "))) %>%
#keep only columns we will use for plotting
select(tip, labels, cor_av, pglsres.cdmass)
# set row names in dataset to match the tip labels in the tree
row.names(sal.phy) <- sal.phy$tip
#check that phylogeny and data match exactly
name.check(caudatatree, sal.phy)
## [1] "OK"
#resort dataset to the order of tree tip labels
sal.phy <- sal.phy[caudatatree$tip.label, ]
#labels
labs <- sal.phy %>% select(tip, labels)
Next, we use ggtree to plot the phylogeny alongside bar plots for absolute cornea diameter and relative conea size.
# Make the phylogeny plot
p <- ggtree(caudatatree) %<+% labs +
geom_tiplab(size = 3, aes(label = labels), fontface = 3) +
xlim_tree(375) +
coord_cartesian(clip = 'off')
# Make a second plot next to the phylogeny (bars for cornea diameter)
p2 <- facet_plot(p, panel="Cornea diameter (mm)", data=sal.phy, geom=geom_segment, aes(x=0, xend=cor_av, y=y, yend=y), color = "#BDC3BE", size = 3)
# Make a third plot next to the first two (bars for relative cornea size)
p3 <- facet_plot(p2, panel='Cornea ~ mass residuals', data=sal.phy, geom=geom_segment, aes(x=0, xend=pglsres.cdmass, y=y, yend=y), color = "#BDC3BE", size=3) +
theme_tree2(legend.position=c(.07, .93)) #add scale bars and move legends
#control widths of panels
gt = ggplot_gtable(ggplot_build(p3))
# gtable_show_layout(gt) # will show you the layout - very handy function
# gt # see plot layout in table format
# gt$layout$l[grep('tree', gt$layout$name)] # you want to find the column specific to panel-2
gt$widths[7] = 0.5*gt$widths[5] # reduce column 7
gt$widths[9] = 0.5*gt$widths[5]
#print figure
plot(gt)
#export figure
pdf("../Outputs/tree-figure.pdf", width = 8, height = 22)
grid::grid.draw(gt) # plot with grid draw
dev.off()
## quartz_off_screen
## 2